% This is the Greenwood et al. (1997) model using the Krusell et al. (2000) CES production function.
% Description:	heterogeneous agent model used in CSC business cylce paper
% filename:		RED1.m 
% files used: 	rbcss, rbclhs, rbcrhs, grad, reorder, solab, vec, qzswitch, hpfilter 
% This version of the model is calibrated as closely as possible to Greenwood et al. (1997).

% Modified to generate the output we need, TvR 04/2012

clear variables;
%clc;

% SET PARAMETER VALUES
global beta;
global deltaeq;
global deltas;
global gamma;
global alfa;
global sk;
global u;
global ws;
global wu;
global req;
global rs;
global y;
global sigma;
global ph;
global mu;
global lambda;
global weight;
global rhoaa;
global rhopp;
global ksshare;
global thetak;
global thetau;
global thetas;
global ces;


% Form of production function; CES = 1, CD = 0
ces = 1;

beta 	= 0.9875;	 %?0.9925   % discount factor
gamma 	= 1;			% inverse of the intertemporal elasticity of substitution
deltas 	= 0.014;		% deprecitaion rate for capital structures
deltaeq	= 0.031;        % depreciation rate for capital equipment

rhoaa 	= 0.9999; %0.95;         % persistence of neutral technology shocks      
rhopp 	= 0.9999; %0.64;         % persistence of investment specific technology shocks

%standard deviation of shocks
%randn('state',4);

sigz  = 0.00765 * sqrt( (1-rhoaa^2) / (1-0.95^2) );
sigq  = 0.0038 * sqrt( (1-rhopp^2) / (1-0.64^2) );  
rhozq = 0; %-.31;
sigzq=rhozq*sigz*sigq;

sk 		= .47;			% skilled share of workforce
u 		= 1-sk;			% unskilled share of workforce
ksshare = 0.13;		    % capital structures share of income
thetak 	= 0.17;         % equipment share in CD case
thetau 	= 0.2887;       % unskilled share in CD case
thetas 	= 1-ksshare-thetak-thetau;  % skilled share in CD case

% ces substitution parameters from Krusell et. al (2000)
sigma 	=  0.401;		% elasticity parameter for u and (ke,s)
ph 		= .685; %-.495;		% elasticity parameter for s and ke

% !! The following parameters must be recalibrated if
% any of the above parameters are changed !! This is done by turning on the
% break on line 119. Then you must iterate on the values below "by hand" until
% you find a unique solution (as described in the paper).

% calibrated parameters 
alfa 	= 0.611; %0.611; 		% share of consumption in utility
weight 	= 0.589; %0.589;        % social planner weight on skilled utility
mu 		= 0.358; %0.413;      	% share parameter on u
lambda 	= .105;  %0.553;      	% share parameter on ke

% VARIABLE NAMEs
disp ('x(1)= capital equpiment')
disp ('x(2)= capital structures')
disp ('x(3)= productivity') 
disp ('x(4)= relative price of capital')
disp ('x(5)= skilled labor supply')
disp ('x(6)= unskilled labor supply')
disp ('x(7)= skilled consumption')
disp ('x(8)= unskilled consumption')

% total number of variables and number of state variables
ns=4;
nd=4;
nvar=nd+ns;

% Find the Steady State

% define a starting point 
x0=[0.5 0.5 0 0 0.5 0.5 0.5 0.00005]';
options(1)=0;
xss=fsolve('rbcss',x0,options);
disp('Steady State'), disp(xss);

keq=xss(1);
ks=xss(2);
z=xss(3);
q=xss(4);
hs=xss(5);
hu=xss(6);
cs=xss(7);
cu=xss(8);

Keq_share=req*keq/y
Ks_share=rs*ks/y;
S_share=sk*hs*ws/y;
U_share=u*hu*wu/y;
wp=ws/wu
ratio=hs/hu 
avg=sk*hs+u*hu
total = Keq_share+Ks_share+U_share+S_share;
K_Share = Ks_share+Keq_share;

%break; 

% Linearize the system around the s.s.
phi=grad('rbcrhs',xss);
psi=grad('rbclhs',xss);

% solution to the dynamic system
% using the generalized complex schur decomposition

[f,p]=solab(psi,phi,ns); 
% f is the decision rules, i.e. the control variables as functions of state variables
% in deviations from the ss
% p gives the law of motion of the state variables as
% functions of themselves - in deviations from the ss
% to get levels, add respectively

disp 'Laws of motion state variables:';
p1=p*xss(1:ns)+(eye(ns)-p)*xss(1:ns)
disp 'Policy rules:';
f1=f*xss(1:ns)+xss(ns+1:nvar)-f*xss(1:ns) 

% simulation - stochastic
s=88; %88; %1000; %96;				% # of periods
t=1:s;
disp 'Number of time periods per simulation:';
s

xs=zeros(nvar,s);	% set dimension of the output matrix
j=1;
J=1; %1000; %300;     	   	    % #of simulations
disp 'Number of simulations:';
J
irf_and_graphs=1;

% Initialize storage vectors 
% Output
yy=zeros(J,1);			
yym1=zeros(J,1);
yym2=zeros(J,1);
yym3=zeros(J,1);
yym4=zeros(J,1);
yym5=zeros(J,1);
yyp1=zeros(J,1);
yyp2=zeros(J,1);
yyp3=zeros(J,1);
yyp4=zeros(J,1);
yyp5=zeros(J,1);
% total capital
yk=zeros(J,1);			
ykm1=zeros(J,1);
ykm2=zeros(J,1);
ykm3=zeros(J,1);
ykm4=zeros(J,1);
ykm5=zeros(J,1);
ykp1=zeros(J,1);
ykp2=zeros(J,1);
ykp3=zeros(J,1);
ykp4=zeros(J,1);
ykp5=zeros(J,1);
% capital equipment
ykeq=zeros(J,1);			
ykeqm1=zeros(J,1);
ykeqm2=zeros(J,1);
ykeqm3=zeros(J,1);
ykeqm4=zeros(J,1);
ykeqm5=zeros(J,1);
ykeqp1=zeros(J,1);
ykeqp2=zeros(J,1);
ykeqp3=zeros(J,1);
ykeqp4=zeros(J,1);
ykeqp5=zeros(J,1);
% capital structures
yks=zeros(J,1);			
yksm1=zeros(J,1);
yksm2=zeros(J,1);
yksm3=zeros(J,1);
yksm4=zeros(J,1);
yksm5=zeros(J,1);
yksp1=zeros(J,1);
yksp2=zeros(J,1);
yksp3=zeros(J,1);
yksp4=zeros(J,1);
yksp5=zeros(J,1);
% consumption
yc=zeros(J,1);			
ycm1=zeros(J,1);
ycm2=zeros(J,1);
ycm3=zeros(J,1);
ycm4=zeros(J,1);
ycm5=zeros(J,1);
ycp1=zeros(J,1);
ycp2=zeros(J,1);
ycp3=zeros(J,1);
ycp4=zeros(J,1);
ycp5=zeros(J,1);
% Skilled consumption
ycs=zeros(J,1);			
ycsm1=zeros(J,1);
ycsm2=zeros(J,1);
ycsm3=zeros(J,1);
ycsm4=zeros(J,1);
ycsm5=zeros(J,1);
ycsp1=zeros(J,1);
ycsp2=zeros(J,1);
ycsp3=zeros(J,1);
ycsp4=zeros(J,1);
ycsp5=zeros(J,1);
% Unskilled consumption
ycu=zeros(J,1);			
ycum1=zeros(J,1);
ycum2=zeros(J,1);
ycum3=zeros(J,1);
ycum4=zeros(J,1);
ycum5=zeros(J,1);
ycup1=zeros(J,1);
ycup2=zeros(J,1);
ycup3=zeros(J,1);
ycup4=zeros(J,1);
ycup5=zeros(J,1);
% Relative consumption
ycscu=zeros(J,1);
ycscum1=zeros(J,1);
ycscum2=zeros(J,1);
ycscum3=zeros(J,1);
ycscum4=zeros(J,1);
ycscum5=zeros(J,1);
ycscup1=zeros(J,1);
ycscup2=zeros(J,1);
ycscup3=zeros(J,1);
ycscup4=zeros(J,1);
ycscup5=zeros(J,1);
% Investment
yi=zeros(J,1);			
yim1=zeros(J,1);
yim2=zeros(J,1);
yim3=zeros(J,1);
yim4=zeros(J,1);
yim5=zeros(J,1);
yip1=zeros(J,1);
yip2=zeros(J,1);
yip3=zeros(J,1);
yip4=zeros(J,1);
yip5=zeros(J,1);
% Productivity
yp=zeros(J,1);			
ypm1=zeros(J,1);
ypm2=zeros(J,1);
ypm3=zeros(J,1);
ypm4=zeros(J,1);
ypm5=zeros(J,1);
ypp1=zeros(J,1);
ypp2=zeros(J,1);
ypp3=zeros(J,1);
ypp4=zeros(J,1);
ypp5=zeros(J,1);
% wage premium
ywp=zeros(J,1);			
ywpm1=zeros(J,1);
ywpm2=zeros(J,1);
ywpm3=zeros(J,1);
ywpm4=zeros(J,1);
ywpm5=zeros(J,1);
ywpp1=zeros(J,1);
ywpp2=zeros(J,1);
ywpp3=zeros(J,1);
ywpp4=zeros(J,1);
ywpp5=zeros(J,1);
% skilled wage
yws=zeros(J,1);			
ywsm1=zeros(J,1);
ywsm2=zeros(J,1);
ywsm3=zeros(J,1);
ywsm4=zeros(J,1);
ywsm5=zeros(J,1);
ywsp1=zeros(J,1);
ywsp2=zeros(J,1);
ywsp3=zeros(J,1);
ywsp4=zeros(J,1);
ywsp5=zeros(J,1);
% unskilled wage
ywu=zeros(J,1);			
ywum1=zeros(J,1);
ywum2=zeros(J,1);
ywum3=zeros(J,1);
ywum4=zeros(J,1);
ywum5=zeros(J,1);
ywup1=zeros(J,1);
ywup2=zeros(J,1);
ywup3=zeros(J,1);
ywup4=zeros(J,1);
ywup5=zeros(J,1);
%relative hours
ysu=zeros(J,1);			
ysum1=zeros(J,1);
ysum2=zeros(J,1);
ysum3=zeros(J,1);
ysum4=zeros(J,1);
ysum5=zeros(J,1);
ysup1=zeros(J,1);
ysup2=zeros(J,1);
ysup3=zeros(J,1);
ysup4=zeros(J,1);
ysup5=zeros(J,1);
% skilled hours
ys=zeros(J,1);			
ysm1=zeros(J,1);
ysm2=zeros(J,1);
ysm3=zeros(J,1);
ysm4=zeros(J,1);
ysm5=zeros(J,1);
ysp1=zeros(J,1);
ysp2=zeros(J,1);
ysp3=zeros(J,1);
ysp4=zeros(J,1);
ysp5=zeros(J,1);
% unskilled hours
yu=zeros(J,1);			
yum1=zeros(J,1);
yum2=zeros(J,1);
yum3=zeros(J,1);
yum4=zeros(J,1);
yum5=zeros(J,1);
yup1=zeros(J,1);
yup2=zeros(J,1);
yup3=zeros(J,1);
yup4=zeros(J,1);
yup5=zeros(J,1);
% total hours
yh=zeros(J,1);			
yhm1=zeros(J,1);
yhm2=zeros(J,1);
yhm3=zeros(J,1);
yhm4=zeros(J,1);
yhm5=zeros(J,1);
yhp1=zeros(J,1);
yhp2=zeros(J,1);
yhp3=zeros(J,1);
yhp4=zeros(J,1);
yhp5=zeros(J,1);
% skilled income
yincs=zeros(J,1);			
yincsm1=zeros(J,1);
yincsm2=zeros(J,1);
yincsm3=zeros(J,1);
yincsm4=zeros(J,1);
yincsm5=zeros(J,1);
yincsp1=zeros(J,1);
yincsp2=zeros(J,1);
yincsp3=zeros(J,1);
yincsp4=zeros(J,1);
yincsp5=zeros(J,1);
% unskilled income
yincu=zeros(J,1);			
yincum1=zeros(J,1);
yincum2=zeros(J,1);
yincum3=zeros(J,1);
yincum4=zeros(J,1);
yincum5=zeros(J,1);
yincup1=zeros(J,1);
yincup2=zeros(J,1);
yincup3=zeros(J,1);
yincup4=zeros(J,1);
yincup5=zeros(J,1);
% skilled investment
yis=zeros(J,1);			
yism1=zeros(J,1);
yism2=zeros(J,1);
yism3=zeros(J,1);
yism4=zeros(J,1);
yism5=zeros(J,1);
yisp1=zeros(J,1);
yisp2=zeros(J,1);
yisp3=zeros(J,1);
yisp4=zeros(J,1);
yisp5=zeros(J,1);
% unskilled investment
yiu=zeros(J,1);			
yium1=zeros(J,1);
yium2=zeros(J,1);
yium3=zeros(J,1);
yium4=zeros(J,1);
yium5=zeros(J,1);
yiup1=zeros(J,1);
yiup2=zeros(J,1);
yiup3=zeros(J,1);
yiup4=zeros(J,1);
yiup5=zeros(J,1);
% capital-skill ratio
yksratio=zeros(J,1);			
yksratiom1=zeros(J,1);
yksratiom2=zeros(J,1);
yksratiom3=zeros(J,1);
yksratiom4=zeros(J,1);
yksratiom5=zeros(J,1);
yksratiop1=zeros(J,1);
yksratiop2=zeros(J,1);
yksratiop3=zeros(J,1);
yksratiop4=zeros(J,1);
yksratiop5=zeros(J,1);
% labor share
ylshare=zeros(J,1);			
ylsharem1=zeros(J,1);
ylsharem2=zeros(J,1);
ylsharem3=zeros(J,1);
ylsharem4=zeros(J,1);
ylsharem5=zeros(J,1);
ylsharep1=zeros(J,1);
ylsharep2=zeros(J,1);
ylsharep3=zeros(J,1);
ylsharep4=zeros(J,1);
ylsharep5=zeros(J,1);
% Storage vectors for variables
y1=zeros(J,1);
c1=zeros(J,1);
cs1=zeros(J,1);
cu1=zeros(J,1);
i1=zeros(J,1);
k1=zeros(J,1);
ks1=zeros(J,1);
keq=zeros(J,1);
h1=zeros(J,1);
hs1=zeros(J,1);
hu1=zeros(J,1);
p1=zeros(J,1);
ws1=zeros(J,1);
wu1=zeros(J,1);
wp1=zeros(J,1);
cr1=zeros(J,1);
hr1=zeros(J,1);
incs1=zeros(J,1);
incu1=zeros(J,1);
is1=zeros(J,1);
iu1=zeros(J,1);
ksratio1=zeros(J,1);
lshare1=zeros(J,1);

% Create array/file for storing simulated data:
VARdata = zeros(s,9,J);
% fid=fopen('VARdata.csv','w');
% fprintf(fid,'Simulated data Lindquist model\n');
% fclose(fid);

while j <= J;
   % draws of shocks
   if rhozq == 0;
      shockz=sigz*randn(1,s+200);
      shockq=sigq*randn(1,s+200);
   else;
      shockz=sigz*randn(1,s+200);
      shockq=sqrt(sigq^2-(sigzq^2/sigz^2))*randn(1,s+200);
      B=sigzq/sigz^2;
      shockq = shockq + mean(shockq) - B*mean(shockz) + B*shockz;
   end;
   shocks(3,:)=shockz; 
   shocks(4,:)=shockq;
   
   % start at the s.s.
   xs(:,1)=xss;
   for i=1:s+200-1;
      xs(1:ns,i+1)=p*xs(1:ns,i)+(eye(ns)-p)*xss(1:ns)+shocks(1:ns,i+1);
   end;
   	xs(ns+1:nvar,:)=f*xs(1:ns,:)+(xss(ns+1:nvar)-f*xss(1:ns))*ones(1,s+200);
      
   keq=	xs(1,:);
   ks	=	xs(2,:);
   k=keq+ks;
   z	=	xs(3,:);
   q	=	xs(4,:);
   hs	=	xs(5,:);
   hu	=	xs(6,:);
   cs	=	xs(7,:);
   cu	=	xs(8,:);
   
   if ces == 1;
      ws=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk*hs).^ph+lambda*keq.^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*...
         (1-ksshare).*(1-mu).*((1-lambda).*(sk*hs).^ph+lambda.*keq.^ph).^((sigma-ph)/ph).*(1-lambda).*(sk.*hs).^(ph-1);
      wu=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk.*hs).^ph+lambda.*keq.^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*...
         (1-ksshare).*mu.*(u.*hu).^(sigma-1);
      y=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk.*hs).^ph+lambda.*keq.^ph).^(sigma/ph)).^((1-ksshare)/sigma);
      rs=exp(z).*ksshare.*ks.^(ksshare-1).*(mu.*(u.*hu).^sigma+(1-mu).*(lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^(sigma/ph)).^((1-ksshare)/sigma);
      req=exp(z).*ks.^ksshare.*((1-ksshare)/sigma).*(mu.*(u.*hu).^sigma+(1-mu).*(lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*(1-mu).*sigma.*...
         (lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^((sigma-ph)/ph).*lambda.*keq.^(ph-1);
   else;
      y = exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^thetau;
      ws = thetas.*exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^(thetas-1).*(u*hu).^thetau;
      wu = thetau.*exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^(thetau-1);
      rs = ksshare.*exp(z).*ks.^(ksshare-1).*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^thetau;
      req = thetak.*exp(z).*ks.^ksshare.*keq.^(thetak-1).*(sk*hs).^thetas.*(u*hu).^thetau;
   end;
   %j
   if floor(j/100)==j/100;
       %j
   end;
   
   % Drop the first 200 observations
   keqtemp=keq(201:s+200);
   keq=keqtemp;
   kstemp=ks(201:s+200);
   ks=kstemp;
   ztemp=z(201:s+200);
   z=ztemp;
   qtemp=q(201:s+200);
   q=qtemp;
   hstemp=hs(201:s+200);
   hs=hstemp;
   hutemp=hu(201:s+200);
   hu=hutemp;
   cstemp=cs(201:s+200);
   cs=cstemp;
   cutemp=cu(201:s+200);
   cu=cutemp;
   c	=	sk*cs+u*cu;
   h	=	sk*hs+u*hu;
   k=keq+ks;
   wutemp=wu(201:s+200);
   wu=wutemp;
   wstemp=ws(201:s+200);
   ws=wstemp;
   ytemp=y(201:s+200);
   y=ytemp;
   rstemp=rs(201:s+200);
   rs=rstemp;
   reqtemp=req(201:s+200);
   req=reqtemp;

   hrsu = u.*hu;
   hrss = sk.*hs;
   relhrs = hrss./hrsu;
   prod=y./h;
   hratio=hu./hs;
   ksratio=keq./hs;
   wp=ws./wu;
   lshare = (sk.*hs.*ws+u.*hu.*wu)./y;
   In= (y-c);

   % Write out the data for the VAR estimates:
   VARdata(:,:,j) = [h' prod' 1./exp(q') wp' wu' ws' relhrs' hrsu' hrss'];
%    fid=fopen('VARdata.csv','a');
%    fprintf(fid,'%s\n','h,prod,q,wp,wu,ws,relhrs,hrsu,hrss');
%    fprintf(fid,'%12.8f,%12.8f,%12.8f,%12.8f,%12.8f,%12.8f,%12.8f,%12.8f,%12.8f\n',[h' prod' exp(q') wp' wu' ws' relhrs' hrsu' hrss']');
%    fclose(fid);
%    %figure(1); plot([h' prod' q' wp' wu' ws' relhrs' hrsu' hrss']);
   
   if irf_and_graphs
       
   % HP-Filter the data
   HP=1;
   if HP==1;
      yt=hpfilter(y,1600);
      kt=hpfilter(k,1600);
      kst=hpfilter(ks,1600);
      keqt=hpfilter(keq,1600);
      zt=hpfilter(z,1600);
      hst=hpfilter(hs,1600);
      hut=hpfilter(hu,1600);
      ht=hpfilter(h,1600);
      hratiot=hpfilter(hratio,1600);
      ct=hpfilter(c,1600);
      Int=hpfilter(In,1600);
      prodt=hpfilter(prod,1600);
      wst=hpfilter(ws,1600);
      wut=hpfilter(wu,1600);
      wpt=hpfilter(wp,1600);
      reqt=hpfilter(req,1600);
      ksratiot=hpfilter(ksratio,1600);
      lsharet=hpfilter(lshare,1600);
      y=y-yt';
      k=k-kt';
      ks=ks-kst';
      keq=keq-keqt';
      z=z-zt';
      hs=hs-hst';
      hu=hu-hut';
      h=h-ht';
      hratio=hratio-hratiot';
      c=c-ct';
      ksratio=ksratio-ksratiot';
      In=In-Int';
      prod=prod-prodt';
      wp=wp-wpt';
      ws=ws-wst';
      wu=wu-wut';
      req=req-reqt';
      lshare=lshare-lsharet';
   end;
   
% calculate percentual fluctuations around trend
y2=y./yt';
In2=In./Int';
k2=k./kt';
keq2=keq./keqt';
ks2=ks./kst';
z2=z./zt';
hs2=hs./hst';
hu2=hu./hut';
h2=h./ht';
hratio2=hratio./hratiot';
ksratio2=ksratio./ksratiot';
c2=c./ct';
In2=In./Int';
prod2=prod./prodt';
ws2=ws./wst';
wu2=wu./wut';
wp2=wp./wpt';
req2=req./reqt';
lshare2=lshare./lsharet';

% calculate and store volatility measures
y1(j,1)=std(y2);
c1(j,1)=std(c2);
i1(j,1)=std(In2);
k1(j,1)=std(k2);
ks1(j,1)=std(ks2);
keq1(j,1)=std(keq2);
h1(j,1)=std(h2);
hs1(j,1)=std(hs2);
hu1(j,1)=std(hu2);
p1(j,1)=std(prod2);
ws1(j,1)=std(ws2);
wu1(j,1)=std(wu2);
wp1(j,1)=std(wp2);
hr1(j,1)=std(hratio2);
ksratio1(j,1)=std(ksratio2);
lshare1(j,1)=std(lshare2);

% corrcoef y:y
   x=corrcoef(y,y);
   yy(j,1)=x(1,2);
   x=corrcoef(y(2:s),y(1:s-1));
   yym1(j,1)=x(1,2);
   x=corrcoef(y(3:s),y(1:s-2));
   yym2(j,1)=x(1,2);
   x=corrcoef(y(4:s),y(1:s-3));
   yym3(j,1)=x(1,2);
   x=corrcoef(y(5:s),y(1:s-4));
   yym4(j,1)=x(1,2);
   x=corrcoef(y(6:s),y(1:s-5));
   yym5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),y(2:s));
   yyp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),y(3:s));
   yyp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),y(4:s));
   yyp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),y(5:s));
   yyp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),y(6:s));
   yyp5(j,1)=x(1,2);
   % corrcoef y:k
   x=corrcoef(y,k);
   yk(j,1)=x(1,2);
   x=corrcoef(y(2:s),k(1:s-1));
   ykm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),k(1:s-2));
   ykm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),k(1:s-3));
   ykm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),k(1:s-4));
   ykm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),k(1:s-5));
   ykm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),k(2:s));
   ykp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),k(3:s));
   ykp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),k(4:s));
   ykp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),k(5:s));
   ykp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),k(6:s));
   ykp5(j,1)=x(1,2);
   % corrcoef y:ks
   x=corrcoef(y,ks);
   yks(j,1)=x(1,2);
   x=corrcoef(y(2:s),ks(1:s-1));
   yksm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),ks(1:s-2));
   yksm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),ks(1:s-3));
   yksm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),ks(1:s-4));
   yksm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),ks(1:s-5));
   yksm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),ks(2:s));
   yksp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),ks(3:s));
   yksp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),ks(4:s));
   yksp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),ks(5:s));
   yksp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),ks(6:s));
   yksp5(j,1)=x(1,2);
   % corrcoef y:keq
   x=corrcoef(y,keq);
   ykeq(j,1)=x(1,2);
   x=corrcoef(y(2:s),keq(1:s-1));
   ykeqm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),keq(1:s-2));
   ykeqm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),keq(1:s-3));
   ykeqm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),keq(1:s-4));
   ykeqm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),keq(1:s-5));
   ykeqm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),keq(2:s));
   ykeqp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),keq(3:s));
   ykeqp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),keq(4:s));
   ykeqp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),keq(5:s));
   ykeqp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),keq(6:s));
   ykeqp5(j,1)=x(1,2);
   % corrcoef y:c
   x=corrcoef(y,c);
   yc(j,1)=x(1,2);
   x=corrcoef(y(2:s),c(1:s-1));
   ycm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),c(1:s-2));
   ycm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),c(1:s-3));
   ycm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),c(1:s-4));
   ycm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),c(1:s-5));
   ycm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),c(2:s));
   ycp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),c(3:s));
   ycp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),c(4:s));
   ycp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),c(5:s));
   ycp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),c(6:s));
   ycp5(j,1)=x(1,2);
   % corrcoef y:invest
   x=corrcoef(y,In);
   yIn(j,1)=x(1,2);
   x=corrcoef(y(2:s),In(1:s-1));
   yInm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),In(1:s-2));
   yInm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),In(1:s-3));
   yInm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),In(1:s-4));
   yInm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),In(1:s-5));
   yInm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),In(2:s));
   yInp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),In(3:s));
   yInp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),In(4:s));
   yInp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),In(5:s));
   yInp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),In(6:s));
   yInp5(j,1)=x(1,2);
   % corrcoef y:z
   x=corrcoef(y,z);
   yp(j,1)=x(1,2);
   x=corrcoef(y(2:s),prod(1:s-1));
   ypm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),prod(1:s-2));
   ypm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),prod(1:s-3));
   ypm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),prod(1:s-4));
   ypm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),prod(1:s-5));
   ypm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),prod(2:s));
   ypp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),prod(3:s));
   ypp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),prod(4:s));
   ypp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),prod(5:s));
   ypp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),prod(6:s));
   ypp5(j,1)=x(1,2);
   % corrcoef y:wp
   x=corrcoef(y,wp);
   ywp(j,1)=x(1,2);
   x=corrcoef(y(2:s),wp(1:s-1));
   ywpm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),wp(1:s-2));
   ywpm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),wp(1:s-3));
   ywpm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),wp(1:s-4));
   ywpm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),wp(1:s-5));
   ywpm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),wp(2:s));
   ywpp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),wp(3:s));
   ywpp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),wp(4:s));
   ywpp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),wp(5:s));
   ywpp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),wp(6:s));
   ywpp5(j,1)=x(1,2);
   % corrcoef y:ws
   x=corrcoef(y,ws);
   yws(j,1)=x(1,2);
   x=corrcoef(y(2:s),ws(1:s-1));
   ywsm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),ws(1:s-2));
   ywsm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),ws(1:s-3));
   ywsm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),ws(1:s-4));
   ywsm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),ws(1:s-5));
   ywsm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),ws(2:s));
   ywsp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),ws(3:s));
   ywsp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),ws(4:s));
   ywsp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),ws(5:s));
   ywsp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),ws(6:s));
   ywsp5(j,1)=x(1,2);
   % corrcoef y:wu
   x=corrcoef(y,wu);
   ywu(j,1)=x(1,2);
   x=corrcoef(y(2:s),wu(1:s-1));
   ywum1(j,1)=x(1,2);
   x=corrcoef(y(3:s),wu(1:s-2));
   ywum2(j,1)=x(1,2);
   x=corrcoef(y(4:s),wu(1:s-3));
   ywum3(j,1)=x(1,2);
   x=corrcoef(y(5:s),wu(1:s-4));
   ywum4(j,1)=x(1,2);
   x=corrcoef(y(6:s),wu(1:s-5));
   ywum5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),wu(2:s));
   ywup1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),wu(3:s));
   ywup2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),wu(4:s));
   ywup3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),wu(5:s));
   ywup4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),wu(6:s));
   ywup5(j,1)=x(1,2);
   % corrcoef y:hs/hu
   x=corrcoef(y,hratio);
   ysu(j,1)=x(1,2);
   x=corrcoef(y(2:s),hratio(1:s-1));
   ysum1(j,1)=x(1,2);
   x=corrcoef(y(3:s),hratio(1:s-2));
   ysum2(j,1)=x(1,2);
   x=corrcoef(y(4:s),hratio(1:s-3));
   ysum3(j,1)=x(1,2);
   x=corrcoef(y(5:s),hratio(1:s-4));
   ysum4(j,1)=x(1,2);
   x=corrcoef(y(6:s),hratio(1:s-5));
   ysum5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),hratio(2:s));
   ysup1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),hratio(3:s));
   ysup2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),hratio(4:s));
   ysup3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),hratio(5:s));
   ysup4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),hratio(6:s));
   ysup5(j,1)=x(1,2);
   % corrcoef y:hs
   x=corrcoef(y,hs);
   yhs(j,1)=x(1,2);
   x=corrcoef(y(2:s),hs(1:s-1));
   yhsm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),hs(1:s-2));
   yhsm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),hs(1:s-3));
   yhsm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),hs(1:s-4));
   yhsm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),hs(1:s-5));
   yhsm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),hs(2:s));
   yhsp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),hs(3:s));
   yhsp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),hs(4:s));
   yhsp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),hs(5:s));
   yhsp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),hs(6:s));
   yhsp5(j,1)=x(1,2);
   % corrcoef y:hu
   x=corrcoef(y,hu);
   yhu(j,1)=x(1,2);
   x=corrcoef(y(2:s),hu(1:s-1));
   yhum1(j,1)=x(1,2);
   x=corrcoef(y(3:s),hu(1:s-2));
   yhum2(j,1)=x(1,2);
   x=corrcoef(y(4:s),hu(1:s-3));
   yhum3(j,1)=x(1,2);
   x=corrcoef(y(5:s),hu(1:s-4));
   yhum4(j,1)=x(1,2);
   x=corrcoef(y(6:s),hu(1:s-5));
   yhum5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),hu(2:s));
   yhup1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),hu(3:s));
   yhup2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),hu(4:s));
   yhup3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),hu(5:s));
   yhup4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),hu(6:s));
   yhup5(j,1)=x(1,2);
   % corrcoef y:h
   x=corrcoef(y,h);
   yh(j,1)=x(1,2);
   x=corrcoef(y(2:s),h(1:s-1));
   yhm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),h(1:s-2));
   yhm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),h(1:s-3));
   yhm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),h(1:s-4));
   yhm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),h(1:s-5));
   yhm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),h(2:s));
   yhp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),h(3:s));
   yhp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),h(4:s));
   yhp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),h(5:s));
   yhp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),h(6:s));
   yhp5(j,1)=x(1,2);
   % corrcoef y:h
   x=corrcoef(y,h);
   yh(j,1)=x(1,2);
   x=corrcoef(y(2:s),h(1:s-1));
   yhm1(j,1)=x(1,2);
   x=corrcoef(y(3:s),h(1:s-2));
   yhm2(j,1)=x(1,2);
   x=corrcoef(y(4:s),h(1:s-3));
   yhm3(j,1)=x(1,2);
   x=corrcoef(y(5:s),h(1:s-4));
   yhm4(j,1)=x(1,2);
   x=corrcoef(y(6:s),h(1:s-5));
   yhm5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),h(2:s));
   yhp1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),h(3:s));
   yhp2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),h(4:s));
   yhp3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),h(5:s));
   yhp4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),h(6:s));
   yhp5(j,1)=x(1,2);
    % corrcoef y:ksratio
   x=corrcoef(y,ksratio);
   yksratio(j,1)=x(1,2);
   x=corrcoef(y(2:s),ksratio(1:s-1));
   yksratiom1(j,1)=x(1,2);
   x=corrcoef(y(3:s),ksratio(1:s-2));
   yksratiom2(j,1)=x(1,2);
   x=corrcoef(y(4:s),ksratio(1:s-3));
   yksratiom3(j,1)=x(1,2);
   x=corrcoef(y(5:s),ksratio(1:s-4));
   yksratiom4(j,1)=x(1,2);
   x=corrcoef(y(6:s),ksratio(1:s-5));
   yksratiom5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),ksratio(2:s));
   yksratiop1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),ksratio(3:s));
   yksratiop2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),ksratio(4:s));
   yksratiop3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),ksratio(5:s));
   yksratiop4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),ksratio(6:s));
   yksratiop5(j,1)=x(1,2);
   % corrcoef y:lshare
   x=corrcoef(y,lshare);
   ylshare(j,1)=x(1,2);
   x=corrcoef(y(2:s),lshare(1:s-1));
   ylsharem1(j,1)=x(1,2);
   x=corrcoef(y(3:s),lshare(1:s-2));
   ylsharem2(j,1)=x(1,2);
   x=corrcoef(y(4:s),lshare(1:s-3));
   ylsharem3(j,1)=x(1,2);
   x=corrcoef(y(5:s),lshare(1:s-4));
   ylsharem4(j,1)=x(1,2);
   x=corrcoef(y(6:s),lshare(1:s-5));
   ylsharem5(j,1)=x(1,2);
   x=corrcoef(y(1:s-1),lshare(2:s));
   ylsharep1(j,1)=x(1,2);
   x=corrcoef(y(1:s-2),lshare(3:s));
   ylsharep2(j,1)=x(1,2);
   x=corrcoef(y(1:s-3),lshare(4:s));
   ylsharep3(j,1)=x(1,2);
   x=corrcoef(y(1:s-4),lshare(5:s));
   ylsharep4(j,1)=x(1,2);
   x=corrcoef(y(1:s-5),lshare(6:s));
   ylsharep5(j,1)=x(1,2);
   
   % if irf_and_graphs
   end
   
   j=j+1;
end;

save VARdata VARdata -v6;

if irf_and_graphs
    
% Volatility output
fid=fopen('Volatility.csv','w');
fprintf(fid,'\n Variable, Volatility (S.D.), S.D. across simulations,');
fprintf(fid,'\n Output ,');
fprintf(fid,'%12.8f,',mean(y1));
fprintf(fid,'%12.8f,',std(y1));
fprintf(fid,'\n capital ,');
fprintf(fid,'%12.8f,',mean(k1));
fprintf(fid,'%12.8f,',std(k1));
fprintf(fid,'\n capital equipment,');
fprintf(fid,'%12.8f,',mean(keq1));
fprintf(fid,'%12.8f,',std(keq1));
fprintf(fid,'\n capital structures,');
fprintf(fid,'%12.8f,',mean(ks1));
fprintf(fid,'%12.8f,',std(ks1));
fprintf(fid,'\n Investment ,');
fprintf(fid,'%12.8f,',mean(i1));
fprintf(fid,'%12.8f,',std(i1));
fprintf(fid,'\n consumption ,');
fprintf(fid,'%12.8f,',mean(c1));
fprintf(fid,'%12.8f,',std(c1));
fprintf(fid,'\n Skilled consumption ,');
fprintf(fid,'%12.8f,',mean(cs1));
fprintf(fid,'%12.8f,',std(cs1));
fprintf(fid,'\n Unskilled consumption ,');
fprintf(fid,'%12.8f,',mean(cu1));
fprintf(fid,'%12.8f,',std(cu1));
fprintf(fid,'\n Hours ,');
fprintf(fid,'%12.8f,',mean(h1));
fprintf(fid,'%12.8f,',std(h1));
fprintf(fid,'\n Unskilled Hours ,');
fprintf(fid,'%12.8f,',mean(hu1));
fprintf(fid,'%12.8f,',std(hu1));
fprintf(fid,'\n Skilled Hours ,');
fprintf(fid,'%12.8f,',mean(hs1));
fprintf(fid,'%12.8f,',std(hs1));
fprintf(fid,'\n Unskilled wages ,');
fprintf(fid,'%12.8f,',mean(wu1));
fprintf(fid,'%12.8f,',std(wu1));
fprintf(fid,'\n Skilled wages ,');
fprintf(fid,'%12.8f,',mean(ws1));
fprintf(fid,'%12.8f,',std(ws1));
fprintf(fid,'\n productivity ,');
fprintf(fid,'%12.8f,',mean(p1));
fprintf(fid,'%12.8f,',std(p1));
fprintf(fid,'\n skill premium ,');
fprintf(fid,'%12.8f,',mean(wp1));
fprintf(fid,'%12.8f,',std(wp1));
fprintf(fid,'\n relative consumption ,');
fprintf(fid,'%12.8f,',mean(cr1));
fprintf(fid,'%12.8f,',std(cr1));
fprintf(fid,'\n relative hours ,');
fprintf(fid,'%12.8f,',mean(hr1));
fprintf(fid,'%12.8f,',std(hr1));
fprintf(fid,'\n skilled income ,');
fprintf(fid,'%12.8f,',mean(incs1));
fprintf(fid,'%12.8f,',std(incs1));
fprintf(fid,'\n unskilled income ,');
fprintf(fid,'%12.8f,',mean(incu1));
fprintf(fid,'%12.8f,',std(incu1));
fprintf(fid,'\n skilled investment ,');
fprintf(fid,'%12.8f,',mean(is1));
fprintf(fid,'%12.8f,',std(is1));
fprintf(fid,'\n unskilled investment ,');
fprintf(fid,'%12.8f,',mean(iu1));
fprintf(fid,'%12.8f,',std(iu1));
fprintf(fid,'\n cap-skill ratio ,');
fprintf(fid,'%12.8f,',mean(ksratio1));
fprintf(fid,'%12.8f,',std(ksratio1));
fprintf(fid,'\n labor share ,');
fprintf(fid,'%12.8f,',mean(lshare1));
fprintf(fid,'%12.8f,',std(lshare1));
fclose(fid);
%type Volatility.txt;

% correlation Output
fid = fopen('Correlations.csv','w');
fprintf(fid,'\n Cross-correlations of Output with:');
A=[yym5 yym4 yym3 yym2 yym1 yy yyp1 yyp2 yyp3 yyp4 yyp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Output:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ykm5 ykm4 ykm3 ykm2 ykm1 yk ykp1 ykp2 ykp3 ykp4 ykp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n capital:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yksm5 yksm4 yksm3 yksm2 yksm1 yks yksp1 yksp2 yksp3 yksp4 yksp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n capital structures:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ykeqm5 ykeqm4 ykeqm3 ykeqm2 ykeqm1 ykeq ykeqp1 ykeqp2 ykeqp3 ykeqp4 ykeqp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n capital equipment:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ycm5 ycm4 ycm3 ycm2 ycm1 yc ycp1 ycp2 ycp3 ycp4 ycp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n consumption:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ycsm5 ycsm4 ycsm3 ycsm2 ycsm1 ycs ycsp1 ycsp2 ycsp3 ycsp4 ycsp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Skilled consumption:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ycum5 ycum4 ycum3 ycum2 ycum1 ycu ycup1 ycup2 ycup3 ycup4 ycup5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Unskilled consumption:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ycscum5 ycscum4 ycscum3 ycscum2 ycscum1 ycscu ycscup1 ycscup2 ycscup3 ycscup4 ycscup5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Relative consumption (cs/cu):,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yInm5 yInm4 yInm3 yInm2 yInm1 yIn yInp1 yInp2 yInp3 yInp4 yInp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Investment:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yhm5 yhm4 yhm3 yhm2 yhm1 yh yhp1 yhp2 yhp3 yhp4 yhp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Hours:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ysum5 ysum4 ysum3 ysum2 ysum1 ysu ysup1 ysup2 ysup3 ysup4 ysup5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Relative Hours (hu/hs):,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yhsm5 yhsm4 yhsm3 yhsm2 yhsm1 yhs yhsp1 yhsp2 yhsp3 yhsp4 yhsp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Skilled Hours:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yhum5 yhum4 yhum3 yhum2 yhum1 yhu yhup1 yhup2 yhup3 yhup4 yhup5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Unskilled Hours:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ywpm5 ywpm4 ywpm3 ywpm2 ywpm1 ywp ywpp1 ywpp2 ywpp3 ywpp4 ywpp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n skill premium (ws/wu):,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ywsm5 ywsm4 ywsm3 ywsm2 ywsm1 yws ywsp1 ywsp2 ywsp3 ywsp4 ywsp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Skilled Wage:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ywum5 ywum4 ywum3 ywum2 ywum1 ywu ywup1 ywup2 ywup3 ywup4 ywup5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n Unskilled Wage:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ypm5 ypm4 ypm3 ypm2 ypm1 yp ypp1 ypp2 ypp3 ypp4 ypp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n productivity:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yincsm5 yincsm4 yincsm3 yincsm2 yincsm1 yincs yincsp1 yincsp2 yincsp3 yincsp4 yincsp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n skilled income:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yincum5 yincum4 yincum3 yincum2 yincum1 yincu yincup1 yincup2 yincup3 yincup4 yincup5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n unskilled income:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yism5 yism4 yism3 yism2 yism1 yis yisp1 yisp2 yisp3 yisp4 yisp5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n skilled investment:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yium5 yium4 yium3 yium2 yium1 yiu yiup1 yiup2 yiup3 yiup4 yiup5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n unskilled investment:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[yksratiom5 yksratiom4 yksratiom3 yksratiom2 yksratiom1 yksratio yksratiop1 yksratiop2 yksratiop3 yksratiop4 yksratiop5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n capital-skill ratio:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
A=[ylsharem5 ylsharem4 ylsharem3 ylsharem2 ylsharem1 ylshare ylsharep1 ylsharep2 ylsharep3 ylsharep4 ylsharep5];
A1=mean(A);
A2=std(A);
fprintf(fid,'\n labor share:,');
fprintf(fid,'%12.8f,',A1);
fprintf(fid,'\n S.D.,');
fprintf(fid,'%12.8f,',A2);
fclose(fid);
%type Correlations.txt;

% break;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% impulse response functions 

% Shocks to Solow residual
t=50;						%periods
out1=zeros(nvar,t);	%set dimension of the output matrix
% draws of shock
shocks=zeros(ns,t);
shocks(3,2)= 0.01;		% 1% one-period technology shock

out1(:,1)=xss;			%start at s.s.
for i=1:t-1
   out1(1:ns,i+1)=p*out1(1:ns,i)+(eye(ns)-p)*xss(1:ns)+shocks(1:ns,i+1);
end
	out1(ns+1:nvar,:)=f*out1(1:ns,:)+(xss(ns+1:nvar)-f*xss(1:ns))*ones(1,t);
   

keq=out1(1,:);
ks=out1(2,:);
z=out1(3,:);
q=out1(4,:);
hs=out1(5,:);
hu=out1(6,:);
cs=out1(7,:);
cu=out1(8,:);
k=ks+keq;
c=sk*cs+u*cu;
h=sk*hs+u*hu;

   if ces == 1;
      ws=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk*hs).^ph+lambda*keq.^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*...
         (1-ksshare).*(1-mu).*((1-lambda).*(sk*hs).^ph+lambda.*keq.^ph).^((sigma-ph)/ph).*(1-lambda).*(sk.*hs).^(ph-1);
      wu=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk.*hs).^ph+lambda.*keq.^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*...
         (1-ksshare).*mu.*(u.*hu).^(sigma-1);
      y=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk.*hs).^ph+lambda.*keq.^ph).^(sigma/ph)).^((1-ksshare)/sigma);
      rs=exp(z).*ksshare.*ks.^(ksshare-1).*(mu.*(u.*hu).^sigma+(1-mu).*(lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^(sigma/ph)).^((1-ksshare)/sigma);
      req=exp(z).*ks.^ksshare.*((1-ksshare)/sigma).*(mu.*(u.*hu).^sigma+(1-mu).*(lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*(1-mu).*sigma.*...
         (lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^((sigma-ph)/ph).*lambda.*keq.^(ph-1);
   else;
      y = exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^thetau;
      ws = thetas.*exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^(thetas-1).*(u*hu).^thetau;
      wu = thetau.*exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^(thetau-1);
      rs = ksshare.*exp(z).*ks.^(ksshare-1).*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^thetau;
      req = thetak.*exp(z).*ks.^ksshare.*keq.^(thetak-1).*(sk*hs).^thetas.*(u*hu).^thetau;
   end;
   
In=y-c;
prod=y./h;
wp=ws./wu;
hratio=hu./hs;
ksratio=keq./hs;
relprice = 1./exp(q);


% transform variables in percentages
ss=ones(nvar,t);
for v=1:nvar ss(v,:)=xss(v); end;
outp=(out1./ss-ones(nvar,t))*100;
y=(y/y(1,1)-ones(1,t))*100;
rs=(rs/rs(1,1)-ones(1,t))*100;
req=(req/req(1,1)-ones(1,t))*100;
ws=(ws/ws(1,1)-ones(1,t))*100;
wu=(wu/wu(1,1)-ones(1,t))*100;
k=(k/k(1,1)-ones(1,t))*100;
ks=(ks/ks(1,1)-ones(1,t))*100;
keq=(keq/keq(1,1)-ones(1,t))*100;
In=(In/In(1,1)-ones(1,t))*100;
h=(h/h(1,1)-ones(1,t))*100;
hs=(hs/hs(1,1)-ones(1,t))*100;
hu=(hu/hu(1,1)-ones(1,t))*100;
wp=(wp/wp(1,1)-ones(1,t))*100;
hratio=(hratio/hratio(1,1)-ones(1,t))*100;
ksratio=(ksratio/ksratio(1,1)-ones(1,t))*100;
c=(c/c(1,1)-ones(1,t))*100;
cs=(cs/cs(1,1)-ones(1,t))*100;
cu=(cu/cu(1,1)-ones(1,t))*100;
prod=(prod/prod(1,1)-ones(1,t))*100;
relprice=(relprice/relprice(1,1)-ones(1,t))*100;

% graphs
tt=linspace(0,t,t);

figure(1);
subplot(3,3,1);
plot(tt,y,'k');
title('output');
ylabel('% deviation');

subplot(3,3,2);
plot(tt,c,'k');
title('consumption');
ylabel('% deviation');

subplot(3,3,3);
plot(tt,In,'k');
title('investment');
ylabel('% deviation');

subplot(3,3,4);
plot(tt,h,'k');
title('hours');
ylabel('% deviation');

subplot(3,3,5);
plot(tt,rs,'k');
title('return on structures');
ylabel('% deviation');

subplot(3,3,6);
plot(tt,req,'k');
title('return on equipment');
ylabel('% deviation');

subplot(3,3,7);
plot(tt,k,'k');
title('total capital');
ylabel('% deviation');

subplot(3,3,8);
plot(tt,ks,'k');
title('capital structures');
ylabel('% deviation');

subplot(3,3,9);
plot(tt,keq,'k');
title('capital equipment');
ylabel('% deviation');

%if ces == 1;
%   print -deps C:\matthew\paper3b\Fig1_CSC.eps;
%else;
%   print -deps C:\matthew\paper3b\Fig1_CD.eps;
%end;

figure(2);
subplot(3,2,1);
plot(tt,hs,'-k',tt,hu,'--k');
legend('skilled   ','unskilled ',1);
ylabel('% deviation');
title('hours');

subplot(3,2,2);
plot(tt,ws,'-k',tt,wu,'--k');
%legend('Skilled','Unskilled',4);
ylabel('% deviation');
title('wages');

subplot(3,2,3);
plot(tt,cs,'-k',tt,cu,'--k');
%legend('Skilled','Unskilled',1);
ylabel('% deviation');
title('consumption');

subplot(3,2,6);
plot(tt,wp,'k');
title('skill premium');
ylabel('% deviation');

subplot(3,2,4);
plot(tt,hratio,'k');
title('relative hours (u/s)');
ylabel('% deviation');

subplot(3,2,5);
plot(tt,ksratio,'k');
title('equipment-skill ratio');
ylabel('% deviation');

%if ces == 1;
%   print -deps C:\matthew\paper3b\Fig3_CSC.eps;
%else;
%   print -deps C:\matthew\paper3b\Fig3_CD.eps;
%end;

figure(5);
subplot(2,4,5);
plot(tt,relprice,'-k');
ylabel('n-shock');
title('price');

subplot(2,4,6);
plot(tt,prod,'-k');
title('productivity');

subplot(2,4,7);
plot(tt,h,'-k');
title('hours');

subplot(2,4,8);
plot(tt,wp,'-k');
title('premium');

% Write out theoretical responses:
theorirf = zeros(t,4,2);  % var, time periods, shock: 1=neutral, 2=i-spec
theorirf(:,:,1) = [relprice',prod',h',wp'];


% Shocks to equipment specific technology
t=50;						%periods
out1=zeros(nvar,t);	%set dimension of the output matrix
% draws of shock
shocks=zeros(ns,t);
shocks(4,2)= 0.01;		% 1% one-period technology shock

out1(:,1)=xss;			%start at s.s.
for i=1:t-1
   out1(1:ns,i+1)=p*out1(1:ns,i)+(eye(ns)-p)*xss(1:ns)+shocks(1:ns,i+1);
end
	out1(ns+1:nvar,:)=f*out1(1:ns,:)+(xss(ns+1:nvar)-f*xss(1:ns))*ones(1,t);
   

keq=out1(1,:);
ks=out1(2,:);
z=out1(3,:);
q=out1(4,:);
hs=out1(5,:);
hu=out1(6,:);
cs=out1(7,:);
cu=out1(8,:);
k=ks+keq;
c=sk*cs+u*cu;
h=sk*hs+u*hu;

   if ces == 1;
      ws=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk*hs).^ph+lambda*keq.^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*...
         (1-ksshare).*(1-mu).*((1-lambda).*(sk*hs).^ph+lambda.*keq.^ph).^((sigma-ph)/ph).*(1-lambda).*(sk.*hs).^(ph-1);
      wu=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk.*hs).^ph+lambda.*keq.^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*...
         (1-ksshare).*mu.*(u.*hu).^(sigma-1);
      y=exp(z).*ks.^ksshare.*(mu.*(u.*hu).^sigma+(1-mu).*((1-lambda).*(sk.*hs).^ph+lambda.*keq.^ph).^(sigma/ph)).^((1-ksshare)/sigma);
      rs=exp(z).*ksshare.*ks.^(ksshare-1).*(mu.*(u.*hu).^sigma+(1-mu).*(lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^(sigma/ph)).^((1-ksshare)/sigma);
      req=exp(z).*ks.^ksshare.*((1-ksshare)/sigma).*(mu.*(u.*hu).^sigma+(1-mu).*(lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^(sigma/ph)).^((1-ksshare-sigma)/sigma).*(1-mu).*sigma.*...
         (lambda.*keq.^ph+(1-lambda).*(sk.*hs).^ph).^((sigma-ph)/ph).*lambda.*keq.^(ph-1);
   else;
      y = exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^thetau;
      ws = thetas.*exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^(thetas-1).*(u*hu).^thetau;
      wu = thetau.*exp(z).*ks.^ksshare.*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^(thetau-1);
      rs = ksshare.*exp(z).*ks.^(ksshare-1).*keq.^thetak.*(sk*hs).^thetas.*(u*hu).^thetau;
      req = thetak.*exp(z).*ks.^ksshare.*keq.^(thetak-1).*(sk*hs).^thetas.*(u*hu).^thetau;
   end;

In=y-c;
prod=y./h;
wp=ws./wu;
hratio=hu./hs;
ksratio=keq./hs;
relprice = 1./exp(q);


% transform variables in percentages
ss=ones(nvar,t);
for v=1:nvar ss(v,:)=xss(v); end;
outp=(out1./ss-ones(nvar,t))*100;
y=(y/y(1,1)-ones(1,t))*100;
rs=(rs/rs(1,1)-ones(1,t))*100;
req=(req/req(1,1)-ones(1,t))*100;
ws=(ws/ws(1,1)-ones(1,t))*100;
wu=(wu/wu(1,1)-ones(1,t))*100;
k=(k/k(1,1)-ones(1,t))*100;
ks=(ks/ks(1,1)-ones(1,t))*100;
keq=(keq/keq(1,1)-ones(1,t))*100;
In=(In/In(1,1)-ones(1,t))*100;
h=(h/h(1,1)-ones(1,t))*100;
hs=(hs/hs(1,1)-ones(1,t))*100;
hu=(hu/hu(1,1)-ones(1,t))*100;
wp=(wp/wp(1,1)-ones(1,t))*100;
hratio=(hratio/hratio(1,1)-ones(1,t))*100;
ksratio=(ksratio/ksratio(1,1)-ones(1,t))*100;
c=(c/c(1,1)-ones(1,t))*100;
cs=(cs/cs(1,1)-ones(1,t))*100;
cu=(cu/cu(1,1)-ones(1,t))*100;
prod=(prod/prod(1,1)-ones(1,t))*100;
relprice=(relprice/relprice(1,1)-ones(1,t))*100;

% graphs
tt=linspace(0,t,t);

figure(3);
subplot(3,3,1);
plot(tt,y,'k');
title('output');
ylabel('% deviation');

subplot(3,3,2);
plot(tt,c,'k');
title('consumption');
ylabel('% deviation');

subplot(3,3,3);
plot(tt,In,'k');
title('investment');
ylabel('% deviation');

subplot(3,3,4);
plot(tt,h,'k');
title('hours');
ylabel('% deviation');

subplot(3,3,5);
plot(tt,rs,'k');
title('return on structures');
ylabel('% deviation');

subplot(3,3,6);
plot(tt,req,'k');
title('return on equipment');
ylabel('% deviation');

subplot(3,3,7);
plot(tt,k,'k');
title('total capital');
ylabel('% deviation');

subplot(3,3,8);
plot(tt,ks,'k');
title('capital structures');
ylabel('% deviation');

subplot(3,3,9);
plot(tt,keq,'k');
title('capital equipment');
ylabel('% deviation');

%if ces == 1;
%   print -deps C:\matthew\paper3b\Fig2_CSC.eps;
%else;
%   print -deps C:\matthew\paper3b\Fig2_CD.eps;
%end;

figure(4);
subplot(3,2,1);
plot(tt,hs,'-k',tt,hu,'--k');
legend('skilled   ','unskilled ',1);
ylabel('% deviation');
title('hours');

subplot(3,2,2);
plot(tt,ws,'-k',tt,wu,'--k');
%legend('Skilled','Unskilled',4);
ylabel('% deviation');
title('wages');

subplot(3,2,3);
plot(tt,cs,'-k',tt,cu,'--k');
%legend('Skilled','Unskilled',1);
ylabel('% deviation');
title('consumption');

subplot(3,2,6);
plot(tt,wp,'k');
title('skill premium');
ylabel('% deviation');

subplot(3,2,4);
plot(tt,hratio,'k');
title('relative hours (u/s)');
ylabel('% deviation');

subplot(3,2,5);
plot(tt,ksratio,'k');
title('equipment-skill ratio');
ylabel('% deviation');

%if ces == 1;
%   print -deps C:\matthew\paper3b\Fig4_CSC.eps;
%else;
%   print -deps C:\matthew\paper3b\Fig4_CD.eps;
%end;

figure(5);
subplot(2,4,1);
plot(tt,relprice,'-k');
ylabel('i-shock');
title('price');

subplot(2,4,2);
plot(tt,prod,'-k');
title('productivity');

subplot(2,4,3);
plot(tt,h,'-k');
title('hours');

subplot(2,4,4);
plot(tt,wp,'-k');
title('premium');

% Write out theoretical responses:
theorirf(:,:,2) = [relprice',prod',h',wp'];
save theorirf theorirf -v6;


% end if irf_and_graphs:
end 
